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SUMMARY 

A key ingredient in the simulation of self-gravitating astrophysical fluid 
dynamical systems is the gravitational potential and its gradient. This paper 
focuses on the development of a mixed method multigrid solver of the Poisson 
equation formulated so that both the potential and the Cartesian components 
of its gradient are self-consistently and accurately generated. The method 
achieves this goal by formulating the problem as a system of four equations 
for the gravitational potential and the three Cartesian components of the 
gradient and solves them using a distributed relaxation technique combined 
with conventional full multigrid V-cyles. The method is described, some 
tests are presented, and the accuracy of the method is assessed. We also 
describe how the method has been incorporated into our three-dimensional 
hydrodynamics code and give an example of an application to the collision of 
two stars. We end with some remarks about the future developments of the 
method and some of the applications in which it will be used in astrophysics. 

1. Introduction 

In recent years a number of astrophysicists [l]-[7] have developed simulation 
tools which build in increasingly realistic physics. The present work grew 
out of an ongoing effort by us to incorporate enough physics and to realize 
that physics with robust algorithms so that we can simulate both existing 
observed phenomena and make reliable predictions which the astronomers 
can utilize in making better observations and interpreting those 
observations. The ubiquitous existence of fluids and gravitation in the 
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universe demands that, if we are to hav e even the most rudimentary 
simulation code, it must incorporate at least interacting fluids and ^ ? r ' 

gravitational physics. In this work, we restrict our attention to the - 
weak-field, Newtonian limit of gr avitation. The hydrodynamics code we 
have created also builds in the effects due to the Special Theory of 
Relativity, so the descriptio n of high sp eed phenomena is included. The 
restriction to weak-field gravity implies that the gr avita tional field is 
determined by the gravitational potential, which must be a solution to 
Poisson’s equation in three dimensions subject to Dirichlet boundary 
conditions at the edges of the computational volume. In the coupled 
hydrodynamic-gravitational system, not only the po ten tial but also its 
gradient is needed. The gradient contributes to the fluid’s acceleration due 
to its self-gravity, inducing the momentum components to change. 

The traditional procedure is to determine the potential by solving the 
Poisson equation with given Dirichlet boundary condition, then construct 
approximations to the components of the gradient via finite differencing the 
potential. However, in simulations of astrophysical gravitating fluids, the 
development of quite complex flows must be anticipated. Examples from 
astrophysics include supernova explosions, gravitational collapse, 
propagation of high-speed jets from active galactic nuclei, star collisions 
and disruptions in dense star clusters, and realistic models of the early 
universe. For most of these simulations, we need to compute the gradients 
of the gravitational potential as accurately as possible, which has motivated 
our development of an alternate appr oach to the gradient computation. 
Here we describe a method which can yield more robust gradients in 
systems that exhibit large variability in space. This is done using a 
distributed relaxation procedur e coupled with full multigri d V- cycles and is 
described in Section 2. In Section 3 we pr esent some tests of the method on 
three-dimensional systems. Section 4 presents our incorporation of it into 
the three-dimensional relativistic hydrodynamics code. Finally, we briefly 
describe an application of the code to the collision of two stars and 
comment on the applications for which the code can be used. 

2. The Mixed Method Algorithm 

The problems we are interested in are three-dimensional, and the results 
that we present in later sections are for such problems. However, in 
presenting the method, we will consider its two-dimensional version to 
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make the description easier to understand and visualize. All components of 
the method, of both the discretization process and the multigrid algorithm, 
have natural three-dimensional analogs. 

a. The Finite Volume Element Discretization Consider the following 
partial differential equation defined on some square domain £1 in H . 

f - V • V</> = / in ft, (i) 

1 <]> = g on dft. 

We let it and v denote the components of the gradient of -<f>: 

(u,vY = — V</> 

Then the partial differential equation may be written in the form of a 
first-order system in ft 

{ u + <f> x = 0 (u equation) 

v + <f> y — 0 (v equation) (2) 

u x + v y = / (p equation), 

with boundary condition 

<f> = g on t?ft. 

Here the labels u,v, and p for the equations are introduced simply for 
convenience. To discretize this system, we follow the Finite Volume 
Element principles developed in [8] . Consider a uniform square mesh ft 
with mesh size h that covers ft. We introduce three sets of control volumes, 
one for each of the three equations in Eq.2. These volumes are shown in 
Fig. 1. We denote by U the set of all volumes U that will be used to 
discretize the u equation in Eq.2. Similarly, we will use the notation V and 
V for the sets of volumes V and P for the v and p equations respectively. 
For our finite element space we consider the lowest order Raviart-Thomas 
elements on the triangulation given by the volumes P: 

u h is linear in x and constant in y on each P € V, 
v h is linear in y and constant in x on each P G P , 

<ft h is constant on each P € P. 

The location of the nodes for each of the unknowns with their indexing is 
also shown in Fig. 1. We can now disretize the equations. We take the u 
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equation in Eq.2 and integrate it over each TJ e U. As an example, let Uf j 
be the volume in U that is centered at the interior u h node We then 
have 


which implies 


/u,j ( u + 0x) dxdy = 0, 

■*“ ^ U i,j U i+lj) M0?+l,i+l — < t > i,j+ 1) = 


Integrating the v equation in Eq.2 over an interior V volume yields a 
similar discrete expression involving nodal values of v h and <f) h . Integrating 
the p equation in Eq.2 over the volume in V centered at the interior <f> h 
node denote this volume by P^; we get 


which implies 


/p t J (u x + v y ) dxdy 


f Pt j fdxdy 




Here, fi tJ is the value of / at the (f> node (i, j) 1 which results from assuming 
that / is (approximated by) a piecewise constant function on V. The only 
remaining part of the discretization involves integrating the u equation in 
Eq.2 over the “half size” U volumes on the left and right boundaries, and 
similarly integrating the v equation in Eq.2 over the “half size” V volumes 
on the lower and upper boundaries. We illustrate this process by 
integrating the u equation in Eq.2 over the volume Ui j that has the 
boundary u h node (1 ,j) as the midpoint of its left edge. We have 


( u + <t>x) dxdy = 0, 

which implies 

^(3 u l ,j + u 2,j) + M02J+1 — = ® 

or 

"h u 2,j) "h d(<f>2 t j+i) = hcfrij+i' 

Note that 4>ij + \ is on the boundary and hence is known. To summarize, the 
discretization has produced for each U volume a discrete version of the u 
equation in Eq.2, for each V volume a discrete version of the v equation in 
Eq.2, and for each P volume a discrete version of the p equation in Eq.2. 

b. The Multigrid Algorithm We assume that the reader is familiar with 
the fundamentals of multigrid methods; good references are [9], [8], [10]. We 
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consider a family of uniform square grids Q h that cover our region ft, where 
h denotes the mesh size. Fig. 2 shows a coarse grid ft? h , with twice the 
mesh size of the grid ft* 1 in Fig. 1. On each grid ft* 1 , we can apply the 
Finite Volume Element discretization process, and we write the discrete set 
of equations that this process generates as 

L h z h = F h j (3) 

where z h = (u h ,v h ,<j> h ) 1 and F h = ( fu h , fv h , f h Y and the unknowns, u h , v h , 
and <j > h , are the nodal values of the corresponding functions on the grid ft . 
Note that the values of 0 at nodes on the boundary are known so they are 
not included in <fi h \ however, as mentioned in the last section, these 
boundary values of d> do appear in the equations generated by integration 
over the U and V volumes near boundaries, resulting in the possibly 
nonzero terms fu h and fv h in Eq.3. In this section, we now define the basic 
components of relaxation, interpolation, and restriction that are necessary 
to implement a multigrid algorithm. 

For the equations on a grid ft* 1 , we use a distributive relaxation process 
similar to that presented in [10]. We can think of relaxation as a three step 
process. First, we sweep over all of the u h nodes, change the value of ttjj so 
that the U equation at (i, j) is satisfied. Second, we perform a similar 
Gauss-Seidel relaxation of all of the V equations. Note that these two steps, 
the U and V relaxation, are independent of each other and could be 
performed in parallel. Finally, we step over the (fY 1 nodes and change the 
value of <f>ij and the values of u h and v h that lie on the edge of the volume 
P*,j, namely and We change these five 

unknowns so that the P equation at (i, j) is satisfied and so that the 
residuals of the U equations at (i,j — 1) and (i — l,j — 1) and of the V 
equations at (i - l t j) and (i - 1, j - 1) are unchanged. To allow 
vectorization, the Gauss-Seidel relaxation performed in each step is done in 
a red/black ordering. 

For defining interpolation operators, we use the same principles as outlined 
in [8] . The Finte Volume Element discretization is based on finite element 
spaces for the variables u h ,v h , and <p k 1 so we can simply use the relationship 
between the finite element spaces on the different grids to define 
interpolation. To define the interpolation operator for <j>, which we denote 
as 7(0)^, we note that <p 2h is constant on the grid 2 h volume P i,j. 
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Referring to Figs. 1 and 2, we thus have the following characterization of 

\ . 7 l . r *-‘ j ;7' 

$,j = $1+1, j = $i,j + 1 = ^i+lj+l ” 0/^j- 


To define the interpolation operator for u, which we denote as I{u)* h , we 
note that u 2h is linear in x and constant in y on the grid 2 h volume P/ j. 
We thus have the following characterization of u h — I [u )\ h u 2h . (See Figs. 1 
and 2) 


u 


= 


- ,.2h 


— U i-1 ,3 — 


U i+l,j-l ~ u i+l,j ~ 
U i,j-1 = u i,j — 


The definition of the interpolation operator for v is similar. 


For defining restriction operators, we again use the same principles as 
outlined in [8]. In the correction scheme multigrid algorithm, which we use 
here, restriction operators are used to transfer right-hand sides and ” — 7 
residuals of equations, not the unknowns themselves. The definitions of the 
restriction operators are based on the relationship between the volumes on 
the various grids. The idea is to lump several of the grid h right-hand sides 
to produce the grid 2h right hand sides. To define the restriction operator 
for the P equation, which we denote as I(P) 2 f, we note that a grid 2 h 
volume P i t j wholly contains four grid h P volumes. We thus have the 
following characterization of f 2h = I(P) 2 h h f h , referring again to Figs. 1 and 
2 : - - 


f I,J = fi,j + f i+l,j + / i,j+l + fi+l,j+l- 

To define the restriction operator for the U equation, which we denote as 
I{U) h > we n °te that a grid 2 h volume U/y in the Interior of fi wholly 
contains two grid h U volumes and half of four others. We thus have the 
following characterization of fu 2h = I(U) 2 ffu h , again referring to Figs. 1 
and 2: 


/«&-• = + Ku-I + + f<i - 1 + /«u + fi&j-J. 


The relationship between U volumes at boundaries is different; for example, 
the grid 2 h U volumes on the left boundary of Q wholly contain two of the 
grid h U volumes and half of two others, yielding Figs. 1 and 2: 


f<j -i = Ki + K,j.i + i/2(/u5j + 
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The definition of the restriction operator for the V equation is done in a 
similar fashion. 

3. Tests of the Mixed Method Algorithm 

A standard approach to Eq.l is to solve a discrete equation based on 
cell-centered finite differences for approximating <fi, then to use simple 
differencing of this approximation to get the components of its gradient. 

We performed some numerical tests to investigate what advantage, in terms 
of accuracy, the mixed method provides over this standard approach. These 
tests were for problems with exact solution 

(j)(x,y,z) = sin(fcix) sin(fc 2 y) sin(fc 3 z) with f 2 = [0, 7r] 3 . By varying ki,k 2 , 
and kz, we were able to see the effect that oscillations in the solution had 
on the accuracy of the methods. Below are results for some of these tests 
on a grid with 32 cells in each direction. 


k , 

k 2 

h 

MIXED METHOD 

STANDARD METHOD 

0err 

(0s) err 

<t>err 

(0x)err 

1 

1 

1 

7.90 E - 4 

8.152? — 4 

1.582? - 3 

8.152? - 4 

1 

16 

16 

1.4715 - 1 

1.502? - 1 

4.592? - 1 

4.722? — 1 | 

16 

1 

1 

1.462? - 1 

3.612? - 0 

4.562? - 1 

3.532? - 0 

16 

16 

16 

1.472? - 1 

3.602? - 0 

rH 

1 

o 

CO 

’3.60 E - 0 


Here, (f) err and {4> x )err are pointwise l 2 norms of the error in <j> and its x 
derivative scaled by the volume term h 3 . These results are indicative of 
results seen for other combinations of fci,fc 2 , and fc 3 . For smooth solutions, 
the methods give nearly identical results. However, for oscillatory solutions, 
the mixed method gives more accurate results, particularly for <j>. 

4. Incorporation of the Mixed Method Solver into 
the Three-Dimensional Hydrodynamics Code 

a.The Physics and the Code The physics included in the present code 
consists of a perfect fluid with an adiabatic equation of state formulated in 
a generally co variant manner. The interval between events in spacetime is 
represented in the present work in the form 

ds 2 = -{a 2 - P^dt 2 + 'Kjidx* + p i dt)(dx j + P j dt). (4) 

The function a is called the lapse and represents the lapse of proper time at 
a given spatial point. The vector field P i is called the shift vector and 



determines how much the spatial coordinates shift from one t — constant 
slice to the next infinitesimally later one. The second rank symmetric 
tensor field 7^ is the metric tensor of the spatial geometry. In the general 
theory of relativity [12], the four- dimensional geometry of spacetime is 
dynamic and the lapse, shift, and three metric are related to the 
kinematical description of the coordinates of the observer and the spatial 
geometry. The fluid energy-momentum tensor must obey a local 
conservation law in order to be consistent with Einstein’s theory. When 
supplemented with the conservation of Baryons, the conservation laws can 
be written in the following form: 

Rest-mass conservation 


1 d 

Q72 dt 


(7 * d ) + 


_L d_ 

cry* dx i 


(7 2 dv l ) - 0 


Internal energy equation 


(5) 


1 d . 1 . 1 d . x ix 

;(7 2 e) + — fai(T ev ) = 


1 Cu ' 

y 2 ut 


ay? ui ay 

j\L 0 m e n t u m equ a t i o n 




ay 


ay 


2 dx l 


( 6 ) 


1 d 


1 d 


dP 


In a 


a 7 


i dt' 


^ s m + = -M +{d+e+ PW)W M 


Si dp* 1 SkSi dl kl 
a dx> 2 W(d + e + PW) dxi 


(7) 


The variables d, e, and 5,, which are used in the code, are defined as 
follows: d = pW , e = peW, and Si = (p + pe + P)v.i . Here d , e, and Si are 
respectively the coordinate mass density, internal energy density, and 
covariant components of the relativistic momentum density. Eq.(5-7) are 
the equations of general relativistic fluid dynamics in a general background 
spacetime. Since the present paper is restricted to the study of phenomena 
with weak gravitational fields, we introduce the following Newtonian 
approximations to the lapse, shift, and three-metric in Cartesian 
coordinates: 


Q ~ l+ 4 > 


( 8 ) 
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F 

in 



(9) 

( 10 ) 


The scalar field <f) is the Newtonian gravitational potential and must satisfy 
the Poisson equation 

VV = 4t rGp, (11) 

in the computational volume and Dirichlet boundary conditions on the 
volume edges. 

With the Newtonian approximation to the geometric variables it then 
follows that the self-gravity of the fluid contributes to the change in the 
momentum density through the term pVa. The value of oc itself enters 
several places in the fluid equations. Thus, a complete characterization of 
the self-gravitating fluid dynamics requires both the lapse and its gradient 
vector. It is these quantities that our mixed method computes in a robust 
manner. Concerning the elements that constitute the hydrodynamics part 
of the code, the methods used may be characterized as explicit finite 
volume schemes. The physical variables d, e, and Si are the fundamental 
quantities. These variables are discretized on a staggered grid system with 
the conventions that scalar variables such as density are stored at zone 
centers, while vector variables are centered on the faces of the zones. The 
biggest challenge is by far to treat the advection of the physical variables as 
accurately as possible. This is especially true for the astrophysical 
applications, since complex flows abound. We want the code to be able to 
detect and track shocks adequately. The advection method implemented in 
the code is based on a monotonic advection algorithm due originally to Van 
Leer [11]. It is robust and tracks shocks reasonably well. The code uses 
artificial viscosity to smooth developing discontinuities over a few zones. 

For this we use an artificial viscosity pressure, which is a combination of 
linear and quadratic functions of the monotonized four-velocity differences. 
The code uses an adiabatic equation of state of the form P = (F — 1 )pe, 
where T is the parameter that characterizes the equation of state and can 
itself be a function of the thermodynamic variables and position. For the 
model stars we discuss here, T is chosen to be a constant. The overall 
structure of a single computational step of the code is described in [7] and 
illustrated as follows: 
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At the end of the computational step the fully updated physical variables 
are available. The Poisson_Solver routine is invoked and it is here that we 
utilize our mixed method solver, which returns (f> and V<j>. 

b. Application to Collision of Stars As a nontrivial application of the 
code, we present a summary of the results of using the mixed method 
Poisson solver in the simulation of the collision of two stars which are 
initially in equilibrium. The initial data were chosen so that the mass 
density and energy density correspond to two equilibrium spheri cal stars. 
WeTiave chosen the n = 1 polytropic equation of state. This equation of 
state has the following functional forms for the initial mass density and 

energy density: d = do and e = eo , where £ = irr/ro and r$ is 

the equilibrium radius of the star. The two model stars were placed with 
their centers displaced in the z = 0 plane. We show here the results of 
simulations in which the radii were chosen equal to 0.26/2 ao j or a nd the, 
central mass density do equal to 6.65/cm 3 . The central temperature of each 
star was chosen to be 4.0e06 K. The simulations shown here were all done 
with a (66) 3 grid. All computations were performed on the Ohio r - 
Supercomputer Center’s Cray YMP8/864. The hydrodynamics part of the 
code has been highly vectorized. 
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Fig. 3a shows the contours for the initial potential and its gradient 
components in the z — 0 plane for a run of an off-center collision. The stars 
were chosen initially to have a relative velocity comparable to the orbital 
velocity. Fig. 3b is a plot of the density contours and velocity field in the 
z=0 plane. Subsequent motion is induced by the combined effects of the 
initial momentum and the self-gravity of the two stars. Because the stars 
attract each other, they develop accelerations toward each other and the 
hydrodynamics that results alters the density and energy distributions. 
Typical simulations were run for at least on the order of the gravitational 
free-fall time. Given the combined interactions of the hydrodynamics with 
self-gravity, we expect disruption of the two stars if the collision is 
sufficiently violent. Figs. 4a, b show respective snapshots of the potential 
contours and gradient and density contours and velocities for late times in 
the off-center collision. 

We conclude from these simulations and others that the mixed method 
Poisson solver produces physically acceptable results when combined with 
the three-dimensional hydrodynamics. This code is currently being used to 
simulate higher resolution runs and other multiple-star systems. We will be 
using the present code to treat the collision of two neutron stars and 
compute its final state and the amount of gravitational radiation emitted 
by such systems. Such computations are of importance because they can 
shed fight on the astrophysics of the mergers of neutron stars as well as 
provide potentially important benchmarks of how much gravitational 
radiation should be expected from such encounters. 
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